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ABSTRACT 



We measure the angular correlation function, w(d), from 0.5 to 30 arcminutes of detected sources in two wide fields of the Herschel Multi-tiered 
Extragalactic Survey (HerMES). Our measurements are consistent with the expected clustering shape from a population of sources that trace the 
dark matter density field, including non-linear clustering at arcminute angular scales arising from multiple sources that occupy the same dark 
^ , matter halos. By making use of the halo model to connect the spatial clustering of sources to the dark matter halo distribution, we estimate source 
"*T* • bias and halo occupation number for dusty sub-mm galaxies at z ~ 2. We find that sub-mm galaxies with 250 ;um flux densities above 30 mjy 
^ ' reside in dark matter halos with mass above (5 ± 4) x IO'^Mq, while (14 ± 8)% of such sources appear as satellites in more massive halos. 
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Cn . 1 . Introduction In the Lockman-SWIRE field we have detected 8154, 4899, 
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and 1680 sources with flux densities above 30 mJy at 250, 350, 
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en ■ The Herschel Multi-tiered Extragalactic Survey (HerMESQ; and 500A/m, respectively, in an area of 218' x 218' (Oliver et al. 

^ ■ Oliver et al. 2010) is the largest project being undertaken by 2010). These counts are supplemented by 3592, 2207, and 1016 

lO Herschel (Pilbratt et al. 2010). It surveys a large set of well- sources detected in the FLS field over an area of 155' x 135', 

O known extra-galactic fields (totaling 70 deg^) at various depths, again down to the same flux density in each of the three bands. 

O primarily using the Spectral and Photometric Imaging Receiver These numbers allow clustering estimates at the same preci- 

^ I (SPIRE) instrument (Griflin et al. 2010). sion level as the first-generation of clustering studies at shorter 

> '. In this Letter, we focus on the clustering of these sources IR wavelengths with source samples from Spitzer data (e.g., 

from 0.5 to 30 arcminute angular scales by making use of source Farrah et al. 2006; Magliocchetti et al. 2007; Waddington et al. 

catalogues in the two widest HerMES fields, Lockman-SWIRE 2007; Brodwin et al. 2008). Instead of simple power-law models, 

j^ and the Spitzer First Look Survey (FLS) observed during the the correlation functions of HerMES sources have high enough 

■ ■ ' Science Demonstration Phase (Oliver et al. 2010). Previous stud- signal-to-noise ratios that we are also able to constrain parame- 

ies on the spatial correlations of sub-mm galaxies was hmited to ters of a halo model (e.g. Cooray & Shefli 2002). 

at most 100 sources, leading either to a limit on the clustering 

amplitude (Blain et al. 2004) or a marginal detection (Scott et 

al. 2006). While the BLAST source catalog was not used for a 2. w{6) measurement 

measurement of the angular correlation function, clustered flue- 

tuations were detected in a power spectrum analysis of all three The angular correlation function, w(6), is a measure of flie prob- 

bands (Viero et al 2009) ability above Poisson fluctuations of finding two galaxies with 

a separation 0, PdD.idD.2 - N[l + w(6)]dQ.[dQ.2, where A^ is 

the surface density of galaxies and dQ.i are solid angles for each 



* Herschel is an ES A space observatory with science instruments pro- galaxy, corresponding to angle 0. The angular coiTelation func- 

vided by European-led Principal Investigator consortia and with impor- tion is of great interest in cosmology as sources are expected to 

tant participation from NASA. trace the underlying dark matter distribution and the clustering 

' hermes.sussex.ac.uk of sources can be related to that of the dark matter halos. 
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Fig. 1. Correlation matrix of the angular correlation function, 

CijI y[Q;C~j, where dj = {w{6i)w(6j)) - <w(0,)><w(0,)) at two 
different angular scales 0, and 6j. Here we show an example 
case at 350 //m for sources in the Lockman-SWIRE field with 
flux densities greater than 30 mJy. 



For clustering measurements we make use of the HerMES 
source catalogues (Oliver et al. 2010), based on maps made 
with calibrated timelines with optimized internal astrometry. The 
maps were produced using the standard SPIRE pipeline after un- 
dergoing calibration and other reduction procedures (Swinyard 
et al. 2010). In the FLS and Lockman-SWIRE fields, a small 
number of individual scans have been removed due to artifacts 
arising from the temperature drift correction. Catalogues were 
generated using the SUSSEXtractor source extractor in HIPE 
version 3.0, using a Gaussian PSF with FWHM of 18.15", 
25.15", and 36.3" at 250, 350, and 500 /im, respectively. Two in- 
dependent maps were also produced by dividing the data in time. 
The sources that do not appear in both sub-maps are flagged as 
spurious sources and removed from the catalogues. The over- 
all astrometry has been adjusted by comparing with radio po- 
sitions. We also apply a Wiener filter optimized for unresolved 
point sources to the maps in the wide fields to remove contam- 
ination from cirrus and then we correct source fluxes through 
simulations. The two fields chosen for this study are largely free 
from cirrus and the clustering measured with two source sam- 
ples with and without the Wiener filter applied agree with each 
other within the errors. Moreover, the source sample used here is 
restricted to sources detected with an overall significance higher 
than 5cr, including confusion noise. 

Once the catalogues are generated, we use the Landy-Szalay 
estimator to measure the correlation function with w{0) - 
[DD(ff) - 2DR(0) + RR(6)]/RR(6), where DD(0) is the num- 
ber of unique pairs of real sources with separation 6, DR{6) is 
the number of unique pairs between the real catalogue and a 
mock sample of sources with random positions, and RR{ff) is 
flie number of unique pairs in the random source catalogues 
(Landy & Szalay 1993). We employ 10^ random mocks with 
10^ sources in each, a larger number of sources than in real 
data to reduce shot-noise in the random pair counts. Since the 
measured w{0) does not automatically satisfy the integral con- 
straint (Infante 1994), we estimate the constant necessary to 
correct w(0) following Adelberger et al. (2005). With a size of 
10 deg^ for each of our two fields, we find a required value of 
~ (2 ± 0.2) X 10"^ between 30' and 40'. This correction is small 
compared to w(0 ~ 30') > 10"^. 
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Fig. 2. Angular correlation function of SPIRE sources in 
Lockman-SWIRE and FLS with flux densities above 30 mJy: (a) 
250 /urn; (b) 350 /im; and (c) 500 fim. The lines are illustrative 
halo models consistent with best-fit results for the occupation 
number (see Table 1), with the dot-dashed lines showing the 2- 
halo term traced by linear clustering and the long-dashed lines 
showing the 1-halo term coming from multiple sources within 
the same halo. The solid lines show the total correlation func- 
tion from our models. 



We calculate the covariance matrix C/j of the correlation 
function involving measurements at two different angular scales 
0i and 0j, using a bootstrap method similar to the one employed 
by Scranton et al. (2002) for measurements of angular cluster- 
ing in SDSS DRl. We also calculate the same covariance ana- 
lytically following the prescription of Eisenstein & Zaldarriaga 
(2001) and find 10% to 20% smaller off-diagonal correlations 
than obtained by bootstrapping the data. We compute the usual 
Poisson errors by taking the square root of the number of pairs, 
finding that the bootstrap variances are a factor of 1.5 to 2 larger 
and are a better representation of errors than the Poisson errors. 
As an example, we show the correlation matrix of 5350 > 30mJy 
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Fig. 3. (a) Angular clustering of 250 /vm sources in Lockman- 
SWIRE field divided into flux densities between 20 and 35 mJy 
and above 35 mJy. (b) Angular clustering for the combined 250 
yum and 350 //m sample again in the Lockman-SWIRE field with 
5350/^250 < 0.85 or > 0.85. 



sources in the Lockman-SWIRE field in Fig. 1, where the corre- 
lation matrix is defined as C,;/ ^CuCjj. 

Since we measure w(6) directly in the real catalogues, it is 
likely to be affected by a variety of effects including source 
blending, flux boosting, and map-making artifacts, among oth- 
ers. To account for all these effects, we compute the transfer 
function necessary to correct w{9) through a set of simulations 
with number counts and source clustering consistent with data 
and in a field similar to Lockman-SWIRE. These input maps are 
then processed by the SPIRE Instrument Simulator (Sibthorpe 
et al. 2009) for an observational program exactly the same as 
HerMES observations of Lockman-SWIRE. The output time- 
ordered data from the simulator are processed as identically as 
the real data. The transfer function is defined as the ratio of the 
correlation function in output catalogues to that of the known 
input used to generate the simulations. The transfer function is 
generally consistent with unity at angular scales of 10' or more, 
but vai-ies 25% (at 250 /im) to 50% (500 /im) at angular scales 
of 3'. We correct the measured w{6) in data based on the average 
of the ratio between input and output w(9) computed from about 
10 simulations. Due to the finite set of simulations we used the 
correction is only known to the level of 20% between 1' and 3' 
angular scales. In comparison, at the same 1' scale, w{9) is un- 
certain at the 30% level at 250 yum and the 50% level at 500 ;um. 
Thus the uncertainty introduced by the error in the transfer func- 
tion is insignificant compared to the overall uncertainties of the 
measured correlation function in data. 
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Fig. 4. Approximate redshift distribution of sources in the 
Lockman-SWIRE field with 5250 > 30 mJy (thick solid line in 
black) and for the two cases based involving colour cuts with 
5350/^250 > 0.85 (magenta dashed line) and 5350/5250 < 0.85 
(orange dot-dashed line) (see text for details). The thin solid 
green line and the two shaded regions in the background show 
example predictions for the S 250 > 30 mJy sample and the two 
colour cuts, respectively, using models from Le Borgne et al. 
(2009; thin green line) and Vahante et al. (2009; shaded regions). 



3. Halo Modeling of Angular Clustering 

In terms of the underlying three-dimensional power spectrum of 
sources as a function of redshift, Ps{k,z), the projected angular 
correlation function is 



w{ff) 



= \ drn\r) f ^P.iK r)Jo(kr6) , 



(1) 



where r(z) is the radial comoving distance, J()(x) is the zeroth 
order Bessel function and n(r) is the radial distribution of sources 
normalized to unity: J drn(r) - 1. 

To model P^ik), we make use of the halo model with both 1- 
and 2-halo terms (Cooray & Sheth 2002). While the 2-halo term 
captures the large-scale clustering with the linear power spec- 
trum scaled by the source bias, the 1-halo term captures the non- 
linear clustering arising from having multiple sources within a 
single dark matter halo. In this case, multiple sources within a 
halo are subdivided to a single source at the halo center and one 
or more satellites. For the central and satellite sources, we de- 
scribe halo occupation numbers as 



<A^ce„(M)) = - 



<A^sat(M)) = - 



1+erf 
l+erf 



log M - log Mmi, 
0"logM 

log M - log 2M„ 

CloeM 



M 



, (2) 



respectively, where erf(x) is the error function, Mmin is the min- 
imum halo mass above which all halos host a central galaxy, 
and the scatter in the relation between galaxy halo mass and 
luminosity is captured by crigg^. We take a fixed value of 0.3 
for (Tiog ^, motivated by modeling of clustering at near-IR wave- 
lengths. Msat is the mass scale at which one satellite galaxy per 
halo is found, in addition to the central galaxy, and a^ is the 
power-law slope of the satellite occupation number with halo 
mass. In addition to these parameters, we also compute the lin- 
ear bias factor of the sources, which should be interpreted as the 
average bias factor of the source sample given the redshift dis- 
tribution (/?;), and the satellite fraction /,, the fraction of sources 
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Table 1. Halo model results using the Lockman-SWIRE w{0) 



Band 


Flux density 


A'gal 


<z> 


log[Mmin/Mo] 


log[M,„/M0] 


Q-i 


<*>c 


/, 


250yum 


5 > 30mJy 


8154 


2 J +04 

9 ri+O'A 
Z.J„7 


12.6+°' 
J 2,9+81 

13.5+^-^ 
13 4+'''-^ 

i2:8!| 


1 3 1 +"■' 


1.3 ±0.4 


2.9 ± 0.4 


0.14 ±0.08 


350//m 


S > 30mJy 


4899 


> 13.1 


< 1.8 


3.2 ±0.5 


<0.20 


SOOjum 


5 > 30mJy 


1680 


2 6+°-^ 


> 13.5 


< 1.6 


3.6 ±0.8 


<0.24 


Combined 


5350/5250 > 0.85 


3333 


2.5 ± 0.4 


> 13.4 


< 1.8 


3.4 ±0.6 


<0.19 


Combined 


5350/^250 < 0.85 


3194 


1 7+"-^ 

^•'-0.6 


> 12.9 


< 1.9 


2.6 ±0.6 


<0.26 



See text below eq. (2) for definitions of M„^i„, Msat, as, (b)., and /j. The redshift range is an approximate estimate based on the colour-colour 
diagram of the source sample through a comparison to isothermal, modified black-body SEDs with a wide range of dust temperatures and emissivity 
parameters (see, Fig. 4 for an example involving 5 250 > 30 mJy and for the two colour cuts). 



in a given sample that appear as satellites in massive dark mat- 
ter halos. This is calculated by taking the ratio of number den- 
sity of satellites to the total number density of sources where 
the number density is calculated through j dM{Ni(M))dn/dM, 
where dn/dM is the halo mass function and index / is for either 
central or satellite source occupation number. 

4. Results & Discussion 

In Fig. 2, we summarize our first set of results related to w{9) 
measurements for each of the three SPIRE bands and for sources 
with S > 30mJy. We show correlation functions measured for 
sources in both Lockman-SWIRE and FLS fields here. We find 
no statistical difference in the correlation functions of sources 
detected in the two fields down to the flux density cut-off of 30 
mJy. To test for evolutionary hints in clustering, in Fig. 3 (a) we 
split the 250 fim source sample to two bins in flux density, while 
in Fig. 3(b) we split the combined sample to two colour bins. 

To model w(6) we need to establish the redshift distribution 
of the source samples. Given the lack of adequate spectroscopic 
redshifts, we make use of sub-mm colours to generate a qual- 
itative redshift distribution (e.g., Hughes et al. 2002). First we 
generate 10^ isothermal SED models using modified black-body 
spectra with a broad range in dust temperature (lOK to 60K) and 
dust emissivity /?. We also include a 15% Gaussian scatter to 
model predictions to account for uncertainties in the observed 
fluxes (Swinyard et al. 2010). We grid the models along the 
redshift direction in the colour-colour plane to several bins and 
simply convert the number of observed data points in each bin 
in the colour-colour plane to a distribution function in redshift 
(see Fig. 4 for an example). While these distributions are gener- 
ally consistent with certain model predictions (e.g., Le Borgne 
et al. 2009; Valiante et al. 2009) and with sub-mm galaxy data 
(e.g.. Chapman et al. 2005), the redshift distribution we recover 
is strongly sensitive to the SEDs used and should only be con- 
sidered as an approximate. 

Our model fitting results are summarized in Table 1 . We as- 
sume WMAP 5-year best-fit ACDM cosmology (Komatsu et al. 
2009). While in Figs. 2 and 3 we only show the errors from 
the variances ( y/Uii), parameter results shown in Table 1 ac- 
count for the covariance matrix when model fitting to measure- 
ments (e.g.. Fig. 1 where we show the correlation matrix). Down 
to the 30 mJy flux density cut, we find average bias factors of 
2.9 ± 0.4 and 3.6 ± 0.8 for 250 and 500 //m sources, respectively. 
Fitting a power-law to all data, the correlation lengths, ro, are 
4.5 ± 0.5 Mpc (250 /im) and 6.3 + 0.7 Mpc (500 jum). While 
250 jum sources are more likely to be found in halos with mass 
(5 + 4) X 10'^ M0, we find that the bright 500 /vm sources in our 
sample occupy halos of (3.1 ± 2.8) x lO'"* Mq. The difference 
is because at a given redshift the 500 yum sources are at a higher 



luminosity. Our modeling allows us to establish that (14 ± 8)% 
of the sources appear as satellites in massive halos than the min- 
imum mass scale. In the case of 350 //m and 500 yum source 
samples, we have failed to accurately determine the parameters 
related to satellite occupation number. As a test on the validity 
of our results to uncertainties in «(z), we also considered two 
extreme possibilities by placing all sources either at z ~ 1 .5 or 
z ~ 3 and found parameters within Icr uncertainties of the esti- 
mates quoted in Table 1 . This is mostly due to the fact that n(z) 
we use for model fitting , with an example shown in Fig. 4, is 
broad with a tail to low redshifts. 

Our measurements show some evidence for non-linear clus- 
tering at arcminute angular scales. Compared to 250 jt/m, an in- 
crease in clustering at arcminute angular scales is less clear at 
350 and 500 //m due to the increase in the beam size. In compar- 
ison, angular power spectra of BLAST fluctuations did not con- 
vincingly reveal a 1-halo term (Viero et al. 2009) and clustering 
was found to be even below the linear term at smallest angular 
scales probed. The increase in arcminute-scale angular cluster- 
ing we see here demonstrates the crucial role played by superior 
angular resolution of SPIRE. Beyond this initial study, future 
work involving understanding the large-scale structure distribu- 
tion of sub-mm galaxies will pursue additional cross-clustering 
studies of sub-mm sources with optical and shorter IR wave- 
lengths, and studies of unresolved fluctuations. On the modeling 
side an approach based on the conditional luminosity function 
could be used to extract additional details on the spatial distribu- 
tion of He rschel sources. 
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